//:Multigrid.cpp
#include "Multigrid.h"


void Restrict::fullweighting(double *f_2h,double *r_h,int m){
    //m :sizeof f_2h
    //std::cout<<"fullweighting_m:"<< m <<std::endl;
    for(int j = 0;j < m;j = j + 1)
        f_2h[j] = r_h[2*j]/4 + r_h[2*j+1]/2 + r_h[2*j+2]/4;
}

void Restrict::injection(double *f_2h,double *r_h,int m){
    //m:sizeof f_2h
    for(int j = 0;j < m;j = j + 1)
        f_2h[j] = r_h[2*j+1];
}

void Restrict::run(double *f_2h,double *r_h,int m){
    if(ISFW) fullweighting(f_2h,r_h,m);
    else injection(f_2h,r_h,m);
}

void Interpolate::linear(double *e_h,double *v_2h,int m){
    //m:sizeof e_h
    e_h[0] =  v_2h[0]/2;
    e_h[m-1] = v_2h[(m-3)/2]/2;
    e_h[m-2] = v_2h[(m-3)/2];
    for(int j = 0;j < (m-3)/2;j = j + 1){
        e_h[2*j+1] = v_2h[j];
        e_h[2*j+2] = (v_2h[j]+v_2h[j+1])/2;
    }
}

void Interpolate::quadratic(double *e_h,double *v_2h,int m){
    //m: sizeof e_h
    e_h[0] = 3*v_2h[0]/4 - v_2h[1]/8;
    e_h[1] = v_2h[0];
    e_h[2] = 3*v_2h[0]/4 + 3*v_2h[1]/8;
    e_h[3] = v_2h[1];
    for(int j = 2;j < (m-3)/2;j = j + 2){
        e_h[2*j] = 3*v_2h[j-1]/8 + 3*v_2h[j]/4 - v_2h[j+1]/8;
        e_h[2*j+1] = v_2h[j];
        e_h[2*j+2] = -v_2h[j-1]/8 + 3*v_2h[j]/4 + 3*v_2h[j+1]/8;
        e_h[2*j+3] = v_2h[j+1];
    }
    e_h[m-3] = 3*v_2h[(m-5)/2]/8 + 3*v_2h[(m-3)/2]/4;
    e_h[m-2] = v_2h[(m-3)/2];
    e_h[m-1] = -v_2h[(m-5)/2]/8 + 3*v_2h[(m-3)/2]/4;
}

void Interpolate::run(double *e_h,double *v_2h,int m){
    if(ISL) linear(e_h,v_2h,m);
    else quadratic(e_h,v_2h,m);
}

void Iteration::weightedJacobi(double *v_h,double *f_h,int m){
    //m:sizeof v_h
    double h = 1.0/(m+1);
    double u[m];
    u[0] = (v_h[1] + h*h*f_h[0]) / 2;
    u[m-1] = (v_h[m-2] + h*h*f_h[m-1]) / 2;
    for(int j = 1;j < m-1;j = j + 1)
        u[j] = (v_h[j-1] + v_h[j+1] + h*h*f_h[j]) / 2;
    for(int j = 0;j < m;j = j + 1)
        v_h[j] = (1-w) * v_h[j] + w * u[j];
}

void Bottomsolver::lapack(double *v_h,double *f_h,int m){
    //m: sizeof v_h
    //std::cout << "begin_solve_m:" << m << std::endl;
    double h = 1.0/(m+1);
    
    double *A_h = new double[m*m];
    double *b = new double[m];

    for(int i = 0;i < m;i = i + 1){
        for(int j = 0;j < m;j = j + 1){
            if(i == j) A_h[i*m+j] = 2/h/h;
            else if(abs(i-j)==1) A_h[i*m+j] = -1/h/h;
            else A_h[i*m+j] = 0;
        }
        b[i] = f_h[i];
    }
    
    int *ipiv = new int[m];
    int info = LAPACKE_dgesv(LAPACK_COL_MAJOR,m,1,A_h,m,ipiv,b,m);
    if(info == 0)
        for(int j =0;j < m;j = j + 1)
            v_h[j] = b[j];
    delete [] A_h;
    delete [] ipiv;
    delete [] b;
    //std::cout << "end_solve_m:" << m << std::endl;
}

double Operator::maxnorm(double *v_h,int m){
    //m:sizeof v_h
    double s = 0;
    for(int j = 0;j < m;j = j + 1)
        if(abs(v_h[j]) > s) s = abs(v_h[j]) ;
    return s;
}

void Operator::residual(double *r_h,double *f_h,double *v_h,int m){
    //m: sizeof v_h
    double h = 1.0/(m+1);
    r_h[0] = f_h[0] + (-2*v_h[0]+v_h[1])/h/h;
    r_h[m-1] = f_h[m-1] + (v_h[m-2]-2*v_h[m-1])/h/h;
    for(int j = 1;j < m-1;j = j + 1)
        r_h[j] = f_h[j] + (v_h[j-1]-2*v_h[j]+v_h[j+1])/h/h;
}

void Operator::correct(double *v_h,double *e_h,int m){
    //m: sizeof v_h
    for(int j = 0;j < m;j = j + 1)
        v_h[j] = v_h[j] + e_h[j];
}
Multigrid::Multigrid(int N){
    n = N;
    l = log(n)/log(2)-2;
    v = new double*[l+1];
    f = new double*[l+1];
    r = new double[n-1];
    e = new double[n-1];
    int m = n;
    for(int j = 0;j < l+1;j = j + 1){
        v[j] = new double[m-1]();
        f[j] = new double[m-1]();
        m = m / 2;
    }
}
void Multigrid::initialize(double *v_h,double *f_h){
    for(int j = 0;j < n-1;j = j + 1){
        v[0][j] = v_h[j];
        f[0][j] = f_h[j];
    }
}

void Multigrid::VCycle(int m,int num){
    //m: sizeof v_h
    //std::cout << "begin_VCycle_m:" << m << std::endl;
    if(num >= l ){
        bottom.lapack(v[num],f[num],m);
        return;
    }
    
    for(int i = 0;i < V1;i = i + 1){
        iteration.weightedJacobi(v[num],f[num],m);
    }

    op.residual(r,f[num],v[num],m);
    op.res.run(f[num+1],r,(m-1)/2);
    for(int j = 0;j < (m-1)/2;j = j + 1)
        v[num+1][j] = 0;
    
    VCycle((m-1)/2,num+1);

    op.inte.run(e,v[num+1],m);
    op.correct(v[num],e,m);
    
    for(int i = 0;i < V2;i = i + 1){
        iteration.weightedJacobi(v[num],f[num],m);
    }
    //std::cout << "end_VCycle_m:" << m << std::endl;
    return;
}

void Multigrid::FullMultigridCycle(int m,int num){
    //m:sizeof v_h
    //std::cout << "Fullmultigridcycle_m:" << m << std::endl;
    if(num >= l){
        VCycle(m,num);return;
    }
    op.res.run(f[num+1],f[num],(m-1)/2);
    
    FullMultigridCycle((m-1)/2,num+1);

    op.inte.run(v[num],v[num+1],m);

    VCycle(m,num);
    return;
}

void Multigrid::solve(double *v0,double *f0,double *u0,int M,double eps,std::ofstream& out){
    //m:sizeof v
    double *v_0 = new double[n-1]();
    double *error = new double[n-1];
    double r0,r1;
    double e0,e1;
    op.res.ISFW = ISFW;
    op.inte.ISL = ISL;
    iteration.w = w;

    out << std::setw(20)<< "iterations number" << std::setw(20) << "residual"
        << std::setw(20) <<"ratio" << std::setw(20) << "error" << std::setw(20) << "ratio" << std::endl;
    op.residual(r,f0,v0,n-1);

    for(int j = 0;j < n-1;j = j + 1)
        error[j] = abs(v0[j]-u0[j]);
    r1 = op.maxnorm(r,n-1);
    e1 = op.maxnorm(error,n-1);
    out << std::setw(20) << "0" << std::setw(20) << std::setiosflags(std::ios::scientific)
        << std::setprecision(5)<< r1 << std::setw(20) << " " << std::setw(20) << e1 << std::setw(20) << " " << std::endl;
    
    for(int i = 0;i < M;i = i + 1){
        if(r1 < eps) break;

        initialize(v_0,r);
        
        if(ISVC)VCycle(n-1,0);
        else FullMultigridCycle(n-1,0);
        
        for(int j = 0;j < n-1;j = j + 1){
            v0[j] = v0[j] + v[0][j];
            error[j] = abs(v0[j] - u0[j]);
        }
        
        op.residual(r,f0,v0,n-1);
        
        e0 = e1;
        e1 = op.maxnorm(error,n-1);
        r0 = r1;
        r1 = op.maxnorm(r,n-1);
        out << std::setw(20) << i+1;
        out << std::setw(20) << std::setiosflags(std::ios::scientific) << std::setprecision(5)
            << r1 << std::resetiosflags(std::ios::scientific) << std::setw(20)
            << std::setiosflags(std::ios::fixed) << std::setprecision(7)
            << r1/r0 << std::resetiosflags(std::ios::fixed) << std::setw(20)
            << std::setiosflags(std::ios::scientific) << std::setprecision(5)
            << e1 << std::resetiosflags(std::ios::scientific) << std::setw(20)
            << std::setiosflags(std::ios::fixed) << std::setprecision(7)
            << e1/e0 << std::resetiosflags(std::ios::fixed) << std::endl;
    }
    delete [] v_0;
    delete [] error;
}

Multigrid::~Multigrid(){
    delete [] v;
    delete [] f;
    delete [] r;
    delete [] e;
}

